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We theoretically and numerically study Bogoliubov-wave turbulence in three-dimensional atomic 
Bose-Einstein condensates with the Gross-Pitaevskii equation, investigating three spectra for the 
macroscopic wave function, the density distribution, and the Bogoliubov-wave distribution. In this 
turbulence, Bogoliubov waves play an important role in the behavior of these spectra, so that we 
call it Bogoliubov-wave turbulence. In a previous study [D. Proment et al., Phys. Rev. A 80, 
051603(R) (2009)], a —3/2 power law in the spectrum for the macroscopic wave function was sug¬ 
gested by using weak wave turbulence theory, but we find that another —7/2 power law appears in 
both theoretical and numerical calculations. Furthermore, we focus on the spectrum for the den¬ 
sity distribution, which can be observed in experiments, discussing the possibility of experimental 
observation. Through these analytical and numerical calculations, we also demonstrate that the 
previously neglected condensate dynamics induced by the Bogoliubov waves is remarkably impor¬ 
tant. 

PACS numbers: 03.75.Kk,05.45.-a 


I. INTRODUCTION 

In nature, complex and chaotic flow universally ap¬ 
pears in various systems from small to large scale, for 
instance, in watercourses, rivers, oceans, and the atmo¬ 
sphere. Such a flow, called turbulence, has attracted con¬ 
siderable attention in various fields because understand¬ 
ing it leads to significant development of fundamental 
and applied sciences such as mathematics, physics, and 
technology. Though turbulence has been studied for a 
long time, however, it is considered to be one of the 
unresolved problems and becomes an important issue in 
modern physics. 

Originally, turbulence has been studied in classical flu¬ 
ids, where it is important to investigate statistical quan¬ 
tities such as the probability density function of the ve¬ 
locity field P, 0 . One of the most significant quantities is 
the kinetic energy spectrum corresponding to a two-point 
correlation function of the velocity field, which is known 
to exhibit the famous Kolmogorov —5/3 power law. This 
Kolmogorov law appears as a result of the constant-flux 
energy cascade in wave number space, and the observa¬ 
tion of this power law is very important in the study of 
turbulence. 

In quantum fluids, turbulence also appears and has 
been investigated for many years M- This turbulence 
is called quantum turbulence (QT), in contrast to clas¬ 
sical turbulence (CT), and was originally developed in 
the field of superfluid helium. QT offers rich phenomena 
such as thermal counterflow Isl, d, Richardson cascade 
I3-[^ Kelvin wave casca de ]l3ni3 |. the bottleneck ef¬ 
fect [ 13 , velocity statistics [Ig, and the reconnection of 
quantized vortices [13: leading to a deep understanding 
of turbulence, and it has been actively studied until to¬ 
day. 


Recently, atomic Bose-Einstein condensates (BECs) 
have been providing a novel stage for turbulence stud¬ 
ies in quantum fluids i i [3- This system has three 
advantages for investigating turbulence in quantum flu¬ 
ids. The first advantage is the experimental visualization 
of the density and vortex distribution [l^ , which gives us 
a large amount of information about turbulence. Such an 
observation is difficult in superfluid helium, though the 
vortices have been visualized recently by use of tracer 
particles El El- The second advantage is the high con¬ 
trollability of the system by optical techniques Ell : which 
enable us to change the shape of the quantum fluid and 
the interaction between particles. This can access new 
problems: the geometrical dependence of turbulence such 
as two-dimensional QT. This kind of study is being ac¬ 
tively investigated. The third advantage is the realization 
of multicomponent BECs such as binary [2(ij | and spinor 
BECs [ 21 ], [ 22 I . The hydrodynamics in this kind of quan¬ 
tum fluid will give us novel physics, and investigations in 
this area have recently begun [23l429l| . Thus, turbulence 
in the atomic BECs can address exotic phenomena not 
found in superfluid helium. 

At present, there is some research on turbulence in 
atomic BECs, and most theoretical and numerical stud¬ 
ies consider turbulence with many quantized vortices. In 
three- and two-dimensional systems, when many quan¬ 
tized vortices are nucleated, the Kolmogorov —5/3 power 
law is numerically confirmed in incompressible kinetic en¬ 
ergy [dol - fsll - Recently, two-dimensional QT has been 
actively investigated, with discussion of features charac¬ 
teristic of two-dimensional systems such as inverse energy 
cascades , vortex clustering 34 - 3^ , r^ative tem¬ 

perature [36| , and decay of vortex number 371. 

In experiments, some groups have succeeded in gen¬ 
erating turbulence in atomic BECs 38l44Ij. where it 
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has been possible to obtain turbulence with many quan¬ 
tized vortices. All these experimental studies focus on 
quantized vortices, investigating the anomalous expan¬ 
sion (^ . vortex dynamics [s^, annihilation of vortices 


41| . and so on. However, the kinetic-energy spectrum 


corresponding to a two-point correlation function of the 
velocity field has not yet been observed, and we cannot 
confirm the Kolmogorov —5/3 power law. This means 
that Kolmogorov turbulence has not been observed. 

Let us elucidate the meaning of Kolmogorov turbu¬ 
lence. In the usual hydrodynamic turbulence (HT), the 
vortices give important structures, whose dynamics gen¬ 
erates the —5/3 power law in the kinetic-energy spectrum 
due to a constant-flux cascade for the kinetic energy. 
However, in other systems, various power laws appear 
as a result of the constant-flux cascade for some quan¬ 
tity such as the wave energy, wave action, and so on. In 
these cases, the common physics for the appearance of 
the power law is the constant-flux cascade. In this pa¬ 
per, we call such turbulence with the constant-flux cas¬ 
cade Kolmogorov turbulence. In the turbulence study, it 
would be very important to observe Kolmogorov turbu¬ 
lence; however, observing the superfluid velocity directly 
in turbulence of atomic BECs is difficult, so we cannot 
confirm whether Kolmogorov turbulence appears as of 
now. 

To shed light on the above problem, we focus on 
weak wave turbulence (WT) with a strong condensate in 
atomic BECs because the density profile of the BEC is 
observable and thus plays an important role in observing 
Kolmogorov turbulence. In this turbulence, Bogoliubov 
waves are significant, so we call it Bogoliubov-wave tur¬ 
bulence in this paper. 

Originally, in classical fluids, WT, which is turbulent 
flow dominated by waves, has been studied 4^. It 
is much different from HT because vortex structures are 
regarded as important in HT [l|, Q . There are weak and 
strong kinds of WT, depending on whether the nonlinear¬ 
ity is weak or strong. In weak WT, a constant-flux cas¬ 
cade of wave energy or wave action generates power-law 
behaviors in the spectrum, and Kolmogorov turbulence 
is known to occur in various wave systems such as those 
involving fluid surfaces (gravity and capillary waves) [4^ , 
magnetic substances (spin waves) , and elastic media 
(acoustic waves) [46l |. 

We expect that it is possible to confirm Kolmogorov 
turbulence in atomic BECs by the observation of the den¬ 
sity profile in Bogoliubov-wave turbulence. In this pa¬ 
per, by applying weak WT theory [H, to the Gross- 
Pitaevskii (GP) equation, we theoretically and numeri¬ 
cally study the power-law behaviors for three spectra of 
a macroscopic wave function, density distribution, and 
Bogoliubov wave distribution, discussing the experimen¬ 
tal possibility for observing Kolmogorov turbulence by 
means of the density spectrum. 

There are some previous works on Bogoliubov-wave 


turbulence in three-dimensional systems T7-4^. In 
previous studies 47, |4^ the —3/2 power law in the 
Bogoliubov-wave energy spectrum was analytically de¬ 
rived, and the equation of the fluctuation was derived 
with the assumption that the condensate function has a 
constant amplitude with a phase rotation induced by the 
chemical potential. Subsequent work has suggested 
the —3/2 power law in the spectrum for the macroscopic 
wave function and a result consistent with this power law 
was numerically obtained. 

We reconsider this Bogoliubov-wave turbulence, ana¬ 
lytically deriving the —7/2 power law in the spectrum for 
the macroscopic wave function and numerically confirm¬ 
ing this power law. We find that the condensate dynam¬ 
ics induced by the fluctuation, which is neglected in the 
previous studies, is important for Bogoliubov-wave tur¬ 
bulence. Furthermore, we focus on the spectrum for the 
density distribution, obtaining the —3/2 power law, and 
discuss the experimental possibility of this power law. 

The article is organized as follows. Section H describes 
the GP equation and the spectra for some quantities. In 
Sec. HI, we apply weak WT theory to the GP equa¬ 
tion, deriving the power laws for the spectra. In Sec. IV, 
we show our numerical result for this Bogoliubov-wave 
turbulence. Section V discusses comparison between the 
previous results and ours, the observation of Kolmogorov 
turbulence, and an inconsistency of our analytical calcu¬ 
lation. Finally, we summarize our study in Sec. VI. 


II. FORMULATION 

We address a one-component BEC in a uniform system 
at zero temperature. This system is well described by the 
macroscopic wave function ip obeying the GP equation 
[l9| given by 

with particle mass m and interaction coefficient g. 

In this paper, we focus on spectra for the macroscopic 
wave function ip, the density distribution p = \ip\'^, and 
the Bogoliubov-wave distribution b. The details of the 
Bogoliubov-wave distribution are defined in Sec. HI B. 
The spectrum for the wave function is defined by 

= ^ E (WfeOP) (2) 

fc —Afc/2<|fci|<fc+A/c/2 

with a resolution of Afc = 27r/L in wave-number space 
and a system size L. The brackets indicate an ensemble 
average. The function ip{k) is the Fourier component 
of the macroscopic wave function calculated by T[ip{r)] 
with J^[-] = J ■ e“*^'’’(iy/L‘^ and spatial dimension d. In 
the same way, we can define the spectra for the density 
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and Bogoliubov-wave distributions as 

= ^ E (3) 

fc —Afc/2<|fci|</i;+Afc/2 

C'b(fc) = ^ E (4) 

fe-Afc/2<|fci|<fc+Afc/2 

with the Fourier components of the density distribution 
p(k) = F[p{r)] and the Bogoliubov-wave distribution 
6(fc). The physical meaning of these spectra is that Cw 
and Cd are correlation functions for the macroscopic wave 
function and density profile, and Cb is the Bogoliubov- 
wave distribution in the wave-number space. 


III. APPLICATION OF WEAK WAVE 
TURBULENCE THEORY 


2 7,2 


ih'tpo—4>{k) = -ih4>{k)—il}d 


h^k 

2m 


+ gpo'k’o 2(j>{k) 


+^*{-k) + 2 ^ ^*{ki)^{k2)6{k + fci - fca) 


kik2 


Y, Hk2)^ik3)S{k - k2 - ks)] ,(9) 




where po is the condensate density IV'oP and (5(-) is the 
Kronecker delta. Substituting Eq. ([8]) into Eq. , we 
obtain the following equation for the fluctuation: 


f) - /j2jL2 _ 


' gpo 


0(fc) + r(-fc) 


+2 Y 0*(fci)0(fe2)<5(fe + fei - fca) 

kik2 


+ Y Hk2)Hk3)Sik - fe2 - fea) 

^2^3 


( 10 ) 


We study weak WT in a uniform BEC with the strong 
condensate, deriving the power exponents of C'w(fc), 
C'd(A:), and C'b(fc). For this purpose, we apply weak WT 
theory to the GP equation O- In this application, there 
is an important points, which is the nonlinear dynamics 
for the condensate. We take the squared terms of the 
fluctuation for the macroscopic wave function, treating 
this nonlinear dynamics. In this case, we can also use 
weak WT theory, obtaining these power exponents. 


Therefore, we obtain Eqs. (jS]) and for the condensate 
and the fluctuation dynamics within the weak nonlinear 
condition. 

For application of weak WT theory, we rewrite Eq. 
(US into the canonical form given by 




dH 

d^*{k )’ 


( 11 ) 


A. Equations of the condensate and the fluctuation 


H = H2 + Ha, 


( 12 ) 


We now derive the equations of the condensate and the 
fluctuation. First, we define the condensate tpo and the 
fluctuation (j) as 

?/; = V’o(l + </'), (5) 

where ipo is defined by 

i’o = ^ j dV. (6) 

From Eqs. m and (III), we can derive 

I (j)dV = 0. (7) 


Equations ([6|) and ([7]) mean that the condensate and the 
fluctuation are the k — 0 and fc ^ 0 Fourier components 
of ipik), respectively. In this paper, we consider a weak 
fluctuation from a strong condensate, assuming the weak 
nonlinear condition \(j)\ <C 1. 

We substitute Eq. ([5]) into the GP equation (IT|), de¬ 
riving equations for ipQ and 4>(k) = J^[(j)\ within second 
order of the fluctuation: 


ih—t/jQ = gpot/jQ 
at 


l + ^(2|^(fci)P + ,^(fci)0(-fei))J,(8) 

fel 


H2 = Y 


h^kf 

2m 


+ gpo 




+ 


^ {4>(ki)(t){-ki) + (j)*{ki)(j)*{-ki)^ 


(13) 


H3 = gpo Y d(ki-k2-ka) 

ki,k2.,k3 

X (^*(fei)0(fe2)^(fc3) + KkiW{k2Wi.ka)) . (14) 

Here, note that the condensate function has a time 
dependence. In previous works (47l - l49l| . the condensate 
function ipQ was assumed to be y/^exp{—ipt/H) with 
condensate density po = Nq/L'^, the fe = 0 particle num¬ 
ber Nq at the initial state, and the chemical potential 
M = gpo- However, we keep the second-order fluctua¬ 
tion terms in Eq. ([8]), so that the time dependence of 
ipo is rather complicated. This term is very important 
for calculating the Bogoliubov-wave distribution, which 
is discussed in Sec. IV B. 
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B. Diagonalization of the Hamitonian and 
Bogoliubov-wave distribution equation 


indicate this, we derive the time development equation 
for pq from Eq. which is given by 


To diagonalize the one-body Hamiltonian H 2 , we use 
the Bogoliubov transformation defined by 

(j){k) = u{k)b{k) + v{k)b*{—k), (15) 

where b{k) is the canonical variable for Bogoliubov waves 
(Bogoliubov-wave distribution) and eb(fc) is the disper¬ 
sion relation for this wave defined by 

eb(fc) = \/eo(fc)(eo(fc) +2gpo) (18) 

with eo(fc) = h^k'^/2m. Applying this transformation to 
the Hamiltonian of Eqs. 0 and 0 leads to 

i?2=^eb(fci)|6(fci)P, (19) 

fcl 

iJa = ^ 5{ki-k2-k-i)V{ki,k2,k-i) 

ki.k2,k3 

X (b*{ki)b{k2)b{k3) + b{ki)b*{k2)b*{k3)j 

+ ^ S{ki+k2 + k3)W{ki,k2,k3) 

ki,k2,k3 

X (b*iki)b*{k 2 )b*ik 3 ) + 6(fci)6(fe2)6(fc3)), (20) 

where V and W are the interaction functions for Bogoli¬ 
ubov waves. These functions are given by 

B(fci, ^2, fca) = gpo(uiU2U3 + V 1 V 2 U 3 

+V 1 U 2 V 3 + V 1 V 2 V 3 + U 1 V 2 U 3 + UlU 2 V^ , (21) 

W{ki,k2,k3) = gpo[uiV2V3+ViU2U^, (22) 

where we use simple expressions such as ui and vi for 
u{ki) and v(ki). 

We must notice that the Bogoliubov coefficients u{k) 
and v{k) have time dependence through the condensate 
density po- Thus, the transformation of Eq. 0 is not 
the canonical transformation. However, as shown in the 
following, the time dependence of these coefficients are 
found to give only small corrections in the weak nonlin¬ 
ear condition, so that the transformation of Eq. 0 ap¬ 
proximately becomes the canonical transformation. To 


^ E [Kki)h{-k3) - b*{k3W{-k,)]. (23) 


Then, the time derivative of the Bogoliubov coefficient 
u{k) becomes 


d du{k)dpo 


(24) 


which means that the value of this derivative is the second 
order of the fluctuation. The same thing is true in the 
time derivative of v{k). As a result, in the weak nonlinear 
condition, the time derivative of Eq. 0 is 


~u{k)^b{k) + v(k)^b*{-k) (25) 


since the third order of the fluctuation is neglected in Eq. 
0 . Thus, the Bogoliubov transformation of Eq. 0 
approximately becomes the canonical transformation. 

Based on the above calculation, the time development 
equation for b{k) becomes 




dH 

db*(k)' 


(26) 


This equation can describe weak nonlinear dynamics for 
Bogoliubov waves. In the following subsection, we apply 
weak WT theory to Eqs. (l2^ and (l26l) . discussing the 
behavior of the spectra for the wave function, density 
distribution, and Bogoliubov-wave distribution. 


C. Kinetic equation for Bogoliubov waves 


In weak WT theory, the kinetic equation for the 
spectrum of the canonical variable b{k) can be derived 
11 El- The spectrum is defined by 


n{k)=[^)\\b{k)\^). (27) 

Applying weak WT theory to Eqs. (l23l) and (l26l) . we 
obtain the following kinetic equation: 

^n(fe) = J (n{k,ki,k2)-'Jl{ki,k2,k) 

—TZ{k 2 ,k,ki)jdkidk 2 , (28) 


TZ{k,ki,k 2 ) = 27r|B(fc, fei, fe 2 )P 
x^d(eb(fc) - eb(fci) - eh{k2))6d{k - ki - k2) 

X {n{ki)n{k 2 ) — n{k 2 )n{k) — n{k)n{ki)) (29) 

with the Dirac delta function ^d(’)- In. this applica¬ 
tion, the contributions from & 162&3 and b'^b^b'^ terms 
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in Eq. ([201 are found to be smaller than those of 
other terms 4^, 0, because these terms accompany 
<5d(eb(fc) + eb(fci) + eb(fc 2 ))- In Bogoliubov waves, the 
argument in this Dirac delta function is never satisfied, 
so that the interaction W{ki,k 2 ,kz) does not appear in 
Eqs. (051) an d (051) . The details of this calculation are 
described in [43|, where a nonlinear canonical transfor¬ 
mation causes the interaction bE(fci, A; 2 , fea) to vanish. 
The transformation does not change the expression for 
y(A;i, fc 2 , fca), but it does change the expression for h{k). 
However, this change gives only a small contribution to 
n{k). 

Note that the application of weak WT theory to Eqs. 
(1551) and (1551) is slightly different from usual cases where 
the dispersion relation Cb and the interaction V do not 
have time dependence. In our case, the condensate den¬ 
sity po depends on time, which makes Cb and V depen¬ 
dent on time. Thus, we cannot use the conventional weak 
WT theory in (^ . but we can derive the kinetic equa¬ 
tion (1281) . The detail of this derivation is described in 
Appendix A. 

We are interested in the large-scale dynamics of Bogoli¬ 
ubov waves because, in experiments, it may be difficult 
to observe small-scale dynamics. Here, to elucidate the 
meaning of large scale, let us note the dispersion relation 
for Bogoliubov waves. From Eq. (ITSl) . it follows that 
the k dependence of the dispersion relation drastically 
changes depending on whether or not the amplitude of 
the wave number is larger than fcb = 2^mgpo/h. In the 
wave-number region smaller than fcb, the dispersion re¬ 
lation becomes phononlike, which is expressed by Hcsk 
with sound velocity Cg = ygpolm^ while in the region 
larger than it behaves like that of a free particle. From 
this property of the dispersion relation, in this paper, 
the large (small) scale means the wave-number region is 
smaller (larger) than fcb- 

Finally, we must comment on the importance of the di¬ 
mension of the kinetic equation in the low-wave-number 
region where eb(fc) has linear dispersion. In this region, 
Bogoliubov-wave turbulence is similar to acoustic turbu¬ 
lence, for which, in a two-dimensional systenn the validity 
of the kinetic equation is controversial Thus, 

in the following, we address three-dimensional systems. 

In the next section, to investigate the behavior of the 
spectra in the low-wave-number region, we calculate an 
expression of the interaction V for the low-wave-number 
limit. 


wave numbers of the same order ^ . This locality is also 
discussed in Appendix D. Thus, the wave-number region 
kj < kh (j = 1, 2,3) is considered in the following. 

In the low-wave-number region, the Bogoliubov coeffi¬ 
cients u{k) and v{k) are approximated as 


u{k) ~ 


1 


! + — + -( 


k 

kh 


4 Vfcb 


\ ^ i / fc y 

) ~4V^J 


(30) 


v{k) ~ 


I kh 


k I f k \'^ I ( k 


kh 4 V fcb / 4 V /cb 


(31) 


We substitute Eqs. dSOD and m into the interaction V 
of Eq. (15T1) . obtaining the following expression: 


l/(fel,fe2,fe3)~ ^ 


kl 


8 Y kik2kz 

1 


(fci - fe - k:i) 
kh 


+^(t! - - kj + k% - k,kj 


+klk3 + k2kl -\- kfks - kikj -|- 12kik2k3) 


■ (32) 


As shown by the two delta functions of Eq. (1551) . 
Bogoliubov waves must satisfy both momentum- and 
energy-conservation laws, which are expressed by fci = 
k 2 + fca and eh{ki) = eb(fc 2 ) -|- eb(fc 3 )- These conserva¬ 
tions lead to the rather simple form of Eq. (l32l) . In the 
low-wave-number region, the dispersion relation becomes 


eh(k) ~ hcsk 





(33) 


As a result, the effective interaction 14 for Bogoliubov 
waves becomes 

Vh{ki,k2,k3) = ^^Vkik2k3. (34) 

2kl'^ 


The derivation of this interaction is described in Ap¬ 
pendix B. This interaction satisfies the scaling law 
14(AA;i, Afe 2 , Afca) = A^/^I4(fci, ^ 2 , fea), which is directly 
related to the power exponent of the spectra. In previous 
studies different canonical variables were used to derive 
the interaction function, which is different from Eq. (IMl) . 
but the scaling law is the same 47. 48|. 


D. Interaction of Bogoliubov waves in the 
low-wave-number region 

We now derive the expression for H(A;i, A; 2 , fca) in the 
low-wave-number region. First, we assume that this in¬ 
teraction is local in the turbulent state, which means that 
the dominant contribution to the interaction comes from 


E. Derivation of power exponents of the spectra 

In a stationary state, we can derive the power expo¬ 
nents of the spectra by using the kinetic equation (|58l) 
with the interaction of Eq. ([34|) . Our derivation makes 
use of the fact that a Bogoliubov-wave energy flux is con¬ 
stant in the wave-number space in the weak WT. 
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We integrate the kinetic equation (1281) with the effec¬ 
tive interaction (l34l) and the dispersion relation (1331) , ob¬ 
taining 


= I{k), 


(35) 


= [kl{k-k,fh{kM) 


-kl{ki-kfh{ki,k) 
-kl{k + kif f2{Kki)]dki, 


(36) 


/i(fc, ki) = 9]i{k — ki) (ji{ki)n{k — ki) — n{k)n{ki) 

—n{k)n{k — A:i)^,(37) 


f 2 {k, ki) = n{ki)n{k) — n{k -|- ki)n{k) 

—n{ki)n{k + ki) (38) 

with the Heaviside step function 0 h(’)- Here we as¬ 
sume the isotropy of n{k) = n(fc), which is reasonable 

in isotropic turbulence. In Appendix C, we show the 
derivation of Eqs. (1551) - (1551) . 

From Eq. (l35l) . we can derive the continuity equation 
of the Bogoliubov-wave energy spectrum given by 

^fb(fc) + ^n(fc)=0, (39) 

where fb(fc) and n(fc) are the Bogoliubov-wave energy 
spectrum and the energy flux defined by 

£h{k) = J eh{k)n{k)dklk 

= A'Kk^ €\j{k)n{k) , (40) 


n(fc) = —47r 


ki 


e-b{ki)I{ki) + n{ki)—e^{ki) dki{.Al) 


In the statistically stationary state, the time deriva¬ 
tives of £h{k) in Eq. (1391) is zero, which means that 
the energy flux n(/c) is independent of the wave number. 
Also, the time derivative of eh{k) is zero because the con¬ 
densate density po is stationary in this state. Assuming 
the power law n{k) oc and applying the transforma¬ 
tion kx = kk to Eq. 611), we obtain 

n(yfc) oc (42) 


which exhibits that the energy flux is constant if the 
power exponent m is 9/2. Thus, we derive the power 
law in the Bogoliubov-wave turbulence: 


(43) 


This power law is found to be the local, so that the 
locality assumption in Sec. IV D is justified, which is 
described in Appendix D. Therefore, the spectrum inte¬ 
grated over the solid angle yields 

C'b(fc) oc (44) 


We can derive the power exponent for Cw by using 
Eq. (I44p . In the low-wave-number region, the relation 
between (/(fc) and b{k) becomes 


(j){k) 


b{k) b*{-k) 

Vk Vk '' 


(45) 


which is obtained from Eq. dig. Then, we assume 
isotropy for 6(fc) and use the approximation to derive 
Eq. (1551) . obtaining 


(|</(A)P)cxi(|&(fc)p). 


(46) 


As a result, in the low-wave-number region, the spectrum 
for the macroscopic wave function exhibits the power law 
behavior given by 


C'w(/c) oc k-^Cbik) 

oc (47) 


In this derivation, we use the Fourier component = 

ipoiSik) + (j){k)). 

For the density spectrum, the following expression for 
^ is useful to derive the power exponent of C^'- 

i>{r) = -\/po + dp{r) 

- V'o (^1 + ^ + *<50(r)^ , (48) 

where ^/>o is written as y^exp(z0o)j and 5p and 59 are 
the density and phase fluctuations, respectively, around 
it. Then, the fluctuations of the wave function and the 
density are obtained from 

</(»’) = + * ^^(’')> (49) 

2/Oo 


5p{r) = po(</>(r)-k(/*(r')), (50) 

from which we can express the Fourier component of 5p 
with the canonical variable b. Performing a similar cal¬ 
culation for Eqs. 65l)”6Zl), we can derive the power law 
in the spectrum for the density distribution: 

C'd(fc) oc fcC'b(fc) 

oc (51) 

We consider that this power law is a candidate for an ex¬ 
perimentally observable quantity to confirm the presence 
of Kolmogorov turbulence. 


n{k) oc k 
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As shown in above, the power exponents of Cw and 
Cd are related to Cb through the behavior of Bogoliubov 
coefficients of Eqs. (HU) and (HZl in the low wave num¬ 
ber region. Thus, the difference of the power exponents 
reflects the property of the collective mode. 

In the next section, by numerically calculating the GP 
equation, we discuss these power laws. 


IV. NUMERICAL RESULTS 



We now present our numerical results for Bogoliubov- 
wave turbulence with the GP equation. One of the main 
results is in good agreement with the power exponents 
of Eqs. (I44L (1471) . and m derived by using weak WT 
theory. 


A. Numerical method 

Our numerical calculation treats a three-dimensional 
BEG without a trapping potential. Time and length are 
normalized by r = fi/gpo and ^ = h/y/2mgpQ, and the 
numerical system size L x L x L is 256^ x 256^ x 256^ 
with spatial resolution dx/^ = 1. Because of this resolu¬ 
tion, the quantized vortex dynamics such as nucleation, 
reconnection, etc. cannot be described correctly since its 



FIG. 1: (Color online) Spatial distribution of the density 
p{x,y,z = 128^) at t/r = (a) 0 , (b) 600, (c) 1200 , and (d) 
2500. The smaller structures are nucleated in turn as time 
passes, providing evidence of the direct energy cascade. 


FIG. 2: (Color online) Time development of the spatially av¬ 
eraged absolute values (j)B. and (j>i for real and imaginary parts 
of the fluctuation (p. These graphs exhibit a system satisfying 
the weak nonlinear condition. 



FIG. 3: (Color online) Time development of the amplitude 
po and phase Oq of the condensate V'o- The macroscopic wave 
function is found not to obey y/poexp{—ipt/h). Particularly, 
the deviation of the phase is large. 


core size is the order of However, we address the tur¬ 
bulence without the vortices, so that it is sufficient to use 
this resolution parameter. In this situation, we numeri¬ 
cally solve the GP equation by using the pseudo-spectral 
method. The time propagation is performed by using a 
fourth-order Runge-Kutta method with time resolution 
dt/r = 5 X 10-3. 

We use an initial state in which energy is injected in the 
large scale to confirm the direct energy cascade solutions 
corresponding to Eq. (Hill . (ITTl) . and m- The initial 
state is prepared by the random numbers as follows: 

+ ( 52 ) 

T[(j)\{k) = 3000(i?i -f iR2)fc'‘?'^exp(-fc2^V0-0016),(53) 

where Ri and R 2 are the random numbers in the range 
[0.5,—0.5). Figure 1(a) shows the density distribution 
in the initial state and demonstrates that the energy is 
injected at the large scale. 

In our calculation, dissipation is phenomenologically 
included by replacing ihd/dt with {i — yfaffc —k\j))hd/dt 
in the Fourier transformed Eq. © [ill l3l| . Here, fcb is 
2^mgpo/h, and 7 is the strength of the dissipation. In 


































our calculation, we use 7 = 0.03. We comment on the 
dissipation region. In this paper, we focus on the low- 
wave-number region, adopting the dissipation working in 
the region larger than fcb- There is arbitrariness in the 
choice of the expression for the dissipation. However, in 
a previous paper , the dissipation was found to work 
in the high-wave-number region at low temperature, so 
that we use this expression for the dissipation. 

In summary, we prepare the unstable initial state for 
the confirmation of the direct cascade, numerically cal¬ 
culating the GP equation with the phenomenological dis¬ 
sipation. In this calculation, we do not force the system, 
so that this is decaying turbulence. 


B. Time development of the condensate and the 
fluctuation 

We numerically confirm that our turbulence satisfies 
the weak nonlinear condition \(j)\ <C 1. Figure 2 shows 
the time dependence of the spatially averaged absolute 
values (/>R and (pi for real and imaginary parts of the 
fluctuation (p, from which the order of the fluctuation \p\ 
is found to be about 0.13. 

Figure 3 shows the time development of the amplitude 
po and phase 9q of the condensate ' 00 ! we find that the 
condensate does not obey y/^exp(—i^t/h). This means 
that the fluctuation term in Eq. (jS]) affects the dynamics 
of the condensate. This effect is very important for cal¬ 
culating the Bogoliubov-wave distribution &(fe), which is 
discussed in Sec. IV C. 


C. Time development of the spectra 

We next numerically calculate the spectra Cb, C^, and 
Cd, finding that the numerical results show good agree¬ 
ment with our analytical results of Eqs. (PI), (pl). and 
dSU. Figure 4 shows the time development of the three 
spectra, demonstrating the direct energy cascade from 
small to large wave numbers. This cascade is also con¬ 
firmed by the density distribution as shown in Figs. IHa)- 
lUd), in which small-scale structure is seen to grow. In 
each spectrum, as time passes, the power-law behaviors 
corresponding to Eq. dm), gZl), and (ED gradually grow, 
and the scaling range finally becomes 0.05 ^ < 0.6. 

The fluctuation in < 0.2 is large, which may be the 
finite-size effect. Particularly, the time development of 
Cb is large. 

We note the calculation of the Bogoliubov-wave distri¬ 
bution Wfeb As pointed out in Sec. IV B, in previous 
works [47H49j, the condensate function is assumed to be 
y/^exp(—ifit/h). When we use this condensate func¬ 
tion and calculate the Bogoliubov-wave distribution, the 
spectrum Cb does not show a —5/2 power law and the 
deviation from this law is much larger. When we use 



0.1 




1 





FIG. 4: (Color online) Time development of the spectra for 
the Bogoliubov-wave distribution Cb, the macroscopic wave 
function Cw, and the density distribution Cd. The spectra 
are averaged over six calculations with different initial noise. 
The direct energy cascade occurs, and the power exponents 
of Eq. dm, dSD, and dm show good agreement with the 
numerical results. 


the condensate function defined by Eq. (ED, however, 
the spectrum Cb shows good agreement with Eq. dm- 
Therefore, it is important to take account of the dynam¬ 
ics of the condensate coupled with the fluctuation. Such 
a condensate has been called a quasicondensate , but 
the condensate dynamics of Eq. (ED has not been treated 
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in previous studies. 

Finally, we comment on the relation between the power 
exponents and the strength of the fluctuation. Due to 
our analytical calculation, if the weak nonlinear condi¬ 
tion is satisfied, the energy flux is small, and the power 
exponents of Eqs. (HU), gzl), and (ED can appear, al¬ 
though we do not numerically confirm the appearance of 
the same power laws in the case of weaker fluctuation 
because it takes a much longer time to do the numer¬ 
ical calculation. On the other hand, if the fluctuation 
is strong, the energy flux is large, and other power ex¬ 
ponents may appear. Such a strong WT is the future 
work. 


V. DISCUSSION 

In this section, we derive the three power exponents 
for Cb, Cw, and Cd. In Sec. V A, we detail the dif¬ 
ference between the power exponents in previous results 
and ours. In Sec. V B, we propose some conditions 
to observe Kolmogorov turbulence in atomic BECs and 
comment on the turbulence theory beyond the weak WT 
one. We consider that our result for Cd may be the key 
to achieve this observation. In Sec. V C, we discuss an 
inconsistency of our analytical calculation in Sec. III. 


A. Comparison between previous results and our 
results 

In a previous study (49l |. the spectrum of the macro¬ 
scopic wave function was suggested to obey a —3/2 power 
law in a three-dimensional system. However, our result 
has another —7/2 power law. 

To discuss the reason for this difference, we review 
the derivation of the —3/2 power law in this previous 
paper [i^. In [i^, Zahkarov and Nazarenko studied 
Bogoliubov-wave turbulence using the two canonical vari¬ 
ables (/? = yfp and n = 2hp0 with ij} = ypexp(i0). To 
study the fluctuation around the strong condensate, they 
took expressions ip = ipQ + ipi and H = Hi, where the 
quantity pq corresponds to the condensate and pi and 
Hi correspond to the fluctuations. By using the following 
canonical transformation for these canonical variables, 
they introduced new canonical variables B{k), diagonal¬ 
izing a one-body Hamiltonian as 

= (54) 

fel 

The canonical variables B{k) were defined by 


with the Fourier components pi{k) = B[pi\ and Ili(fe) = 
Here, note that is different from H 2 of Eq. 
(0. Then, they derived the —3/2 power law in the 
Bogoliubov-energy spectrum defined by 

= ^ eb(fc)(|S(fc)H. (57) 

fc—Afc/2<|fci|<fc+A/c/2 

The previous paper (49l | makes use of this power law to 
derive the power exponent of the spectrum for the macro¬ 
scopic wave function. In this previous paper Ssik) is 
estimated as 

^ E eUfc)(|^(fe)p), (58) 

k-Ak/2<\ki\<k+Ak/2 

with e(,(fc) = gpo + eh{k). Then, in the low-wave-number 
region, e(,(fc) ~ gpo is approximately satisfied, so that 
they derived the —3/2 power law in the spectrum Cw for 
the macroscopic wave function. 

However, owing to the weak fluctuation, the macro¬ 
scopic wave function tp can be expanded as 

By applying the Fourier transformation to Eq. dSi) , we 
obtain 

'4>ik)^Po6{k) + pi{k)+i^l^^. (60) 

ZrupQ 

As a result, in the low-wave-number region, the approx¬ 
imation used in weak WT theory gives 

{\m\^)^-^mkr), (61) 

which obviously leads to the —IjZ power law in Cw be¬ 
cause £h{k) oc k{\B[k)\^)k"^ oc Thus, the deriva¬ 

tion of the —3/2 power law for Cw in the previous work 
[ 4 ^ seems not to be correct. 

B. Experimental possibility of Kolmogorov 
turbulence in atomic BECs 

In atomic BEC experiments, observations can be made 
of the density distribution, which accesses the spectrum 
Cd. This can be the key to confirming the presence 
of Kolmogorov turbulence; if we can observe the —3/2 
power law in this spectrum, then Kolmogorov turbulence 
is present. 

We discuss three points related to observing the —3/2 
power law in Cd. The first point is that one must take a 
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wide scaling region; that is, one must prepare an atomic 
BEC that is much larger than the coherence length. The 
—3/2 power law appears in the low-wave-number region 
k < k\j, but, in the trapped system, there is a limit on 
the system size of the order of the Thomas-Fermi radius 
i?TF- This scale can be the lower limit of the scaling 
region for the —3/2 power law. Thus, the radius i?TF 
should be at least a factor of 10 larger than 27r/fcb. 

The second point concerns how to excite the sys¬ 
tem weakly in the low-wave-number region. The —3/2 
power law is valid under the weak nonlinearity condi¬ 
tion 1^1 «C 1, so that fluctuations of both density and 
phase should be very weak. In our calculation, the den¬ 
sity fluctuation is about 5%, as shown in Fig. 1. However, 
this may be too weak to observe. We must investigate 
the spectrum dependence on the strength of the fluctua¬ 
tion, for instance, the condition at which the prediction 
of weak WT theory is broken; this, however, is a subject 
of future work. Also, it is necessary to excite a large- 
scale modulation because the —3/2 power law results 
from the direct energy cascade. Thus, candidates for 
the excitation are considered to be repulsive or attractive 
Gaussian potentials moving slower than the sound veloc¬ 
ity or a shaking trapping potential with fixed frequency. 
However, the excitation mode induced by these methods 
should depend on size, strength, and the velocity of the 
Gaussian potential or frequency and displacement of the 
shaking, etc. Understanding the excitation methods is 
essential for observing Kolmogorov turbulence in atomic 
BEGs and should be numerically investigated hereafter. 

The third point is that one must take surface effects 
into consideration. In a three-dimensional system, we 
can obtain the density distribution integrated along the 
incident direction of the probe light, so that this distri¬ 
bution contains information of the surface profile as well 
as the bulk. Since our result should be applicable to the 
bulk of the trapped BEC, the surface effect may alter 
the power law of Eq. (I5II1 . If the system size is much 
larger than the damping length of the surface wave into 
the bulk 0 , this power law is expected to be relevant. 

Finally, we comment on turbulence theory beyond 
weak WT theory. Recently, Yoshida and Arimitsu have 
studied strong turbulence with the spectral closure ap¬ 
proximation, deriving the some power exponents (54l |. 
These power exponents have not been confirmed numeri¬ 
cally. These power laws may be another key to confirma¬ 
tion of Kolomogorov turbulence in atomic BECs because, 
in this theory, the fluctuation is large and it is easy to 
observe. 


C. An inconsistency of our analytical calculation 

In Sec. HI and Appendix A, we apply weak WT theory 
to Eqs. (1^ and (1^ . obtaining Eq. (|^ . However, this 
derivation has an inconsistency for the perturbative ex¬ 


pansion in Eqs. ([8]) and (HOI) . In this expansion, we keep 
the second order of the fluctuation in both equations, 
which is not consistent because the right-hand sides of 
Eqs. (IH]) and (fTOll are the first and second order of small¬ 
ness, respectively. To do a reliable expansion, we must 
take the third order of 0 in Eq. (nni), but this calculation 
is very difficult. 

Our numerical calculation exhibits that the condensate 
dynamics is different from y^exp(—so that we 
must take the nonlinear dynamics for In this paper, 
although our calculation contains the inconsistency for 
the perturbative expansion, we treat the coupled non¬ 
linear dynamics for the condensate and the fluctuation. 
Thus, our theory for this coupled dynamics is insuffi¬ 
cient. The construction of the improved theory is the 
future work. In spite of this inconsistency, the analyt¬ 
ical results exhibit good agreement with the numerical 
results. 

As another possibility, our numerical calculation seems 
not to be substantially satisfied with the weak nonlinear 
condition. However, if we use the weaker fluctuation, the 
computational time for the formation of the power law is 
much longer, which cannot be confirmed at present. 


VI. CONCLUSION 

We have theoretically and numerically studied 
Bogoliubov-wave turbulence in a uniform BEC at zero 
temperature by using the GP equation. Our application 
of weak WT theory to the GP equation leads to the three 
power exponents —5/2, —7/2, and —3/2, corresponding 
to the spectra for the Bogoliubov-wave distribution Cb, 
the macroscopic wave function Cw, and the density distri¬ 
bution Cd. By numerically calculating the GP equation, 
we have confirmed the good agreement with these power 
exponents. In our analytical and numerical calculation, 
we found that taking account of the condensate dynam¬ 
ics coupled with the fluctuation, an aspect not treated in 
previous works (UTl - ld^ . was important. 
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Appendix A: Application of weak WT theory 


We apply weak WT theory in [43| to Eqs. (l23l) and 
(l26)l . deriving Eq. (|28l) . 

In this theory, the separation of time scale is impor¬ 
tant. There are two time scales, which are a linear time 
scale Tl ^ 27rh/eb and a nonlinear one Tnl- In weak 
WT, the nonlinear term is very weak, which leads to 
the long nonlinear time scale Tnl ^ Tl. In weak WT 
theory, we make use of this time-scale separation, fil¬ 
tering out the fast linear oscillation. For this purpose, 
let us introduce auxiliary intermediate time T satisfying 
Tnl ^ T ^ Tl,. In the following, we solve Eqs. (l23l) 
and (1261) with a formal perturbative expansion, finding 
the canonical variables 6(fc) at t = T. 

Before the perturbative expansion, we rewrite Eqs. 
(USD and (PSI) with interaction variables c(fc) defined by 


b(k) = ec(fe)e-*“('=), 


(Al) 


I f* 

a{k) = -J eb{k,ti)dti, (A2) 

where we introduce the small parameter e in order to 
elucidate the power of smallness of each term in the per¬ 
turbative expansion. Then, Eqs. (l23l) and (1261) become 


zli^c(fe) = e ^ ^S{k - ki - k 2 )V{k,ki,k 2 ) 

ki ,/C 2 

X Ai(fc, ki,k2)c{ki)c{k2) 
+26{ki — k2 — k)V{ki,k2,k) 


xA2{ki,k2,k)c{ki)c*{k2) , 


(A3) 


^ 29pI . 


(A4) 


with Ai(fc, fei, A; 2 ) = exp[z(a(fc) — a(fci) — a(k 2 ))] and 
A 2 {ki,k 2 ,k) = exp[—i(a(fci) — a{k 2 ) — a{k))]. We do 
not show the expression of A because it is not important. 
In this calculation, as noted in Sec. Ill C, the terms with 
W are neglected. 

We calculate the canonical variable c(fc) and the con¬ 
densate density po at t = T up to the second order of e 
by using the perturbative expansion: 


c(k) = c^^\k) + ec^^\k) + e‘^c^‘^\k) -!-•••, (A5) 


^ AO) I I ^2 (2) I f Afi) 

Po — Po + epo + £ Po + • • • ■ (, Ab j 

From Eqs. (TAil) and (IMl) . the expansion of po is given 

by 

po(T) = po(0) + e^p^^\T) + --: (A7) 

The function T(po) such as V, Ai, and A 2 has po as a 
variable, being expanded in the following: 

FipoiT)) F(po(0)) + (A8) 


( 2 ) 

Thus, the second order of the condensate density Pp 
can be neglected in Eq. (IMI) because such terms lead 
to the third order of e in Eq. (IA3I) . Therefore, the time 
dependence of po vanishes in the application of weak WT 
theory, so that we can use the conventional derivation of 
the kinetic equation in [i^, obtaining Eq. (l28l) . 


Appendix B: Derivation of Eq. (1341) 

The derivation of Eq. (1341) is described. In the low- 
wave-number region, the momentum and energy reso¬ 
nant conditions of the Bogoliubov waves give the par¬ 
ticular class of the wave number, where the wave vec¬ 
tors are nearly parallel to each other because of the lin¬ 
ear dispersion [4^ . Then, the resonant condition 

eb(fci) = eb(fc 2 ) + eb(|fci - ^ 2 !) becomes 

ki-k2 + ^{kf - kl) - \ki - k2\ + ^\ki - fe2|^Pl) 


\ki-k2\^A/{ki-k2Y+kik2e'^, (B2) 

where the angle 9 formed by the two vectors fei and k 2 
is small. Taking the second order of the small quantities, 
we obtain the angle 9 ~ V^iki — k 2 )/k\y [i^, which gives 


|fci-fe 2 | 2 ^ (fci - fc 2 )y 1 -k 

^ {ki-k2){l + ^). (B3) 

Here we use the condition ki > ^2 because of the resonant 
condition. Thus, the effective interaction between the 
Bogoliubov waves can be calculated by V (fci, fc 2 , fci — ^ 2 ) 
of Eq. (l32l) . As a result, we obtain 


Eb(fci,fe 2 ,fel - ^ 2 ) 


■\/fcifc2(fei - k2). (B4) 

2k^'^ 


Appendix C: Derivation of Eq. (1351) 

We show the derivation of Eq. (1351) from Eq. (l28l) . We 
assume the isotropy of n(fc) in the weak WT, integrating 
the k 2 part and taking the leading terms, which leads to 

^n(A)=J(fc), (Cl) 
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2(fc _ k^)S 40 l - ^Dfik, kuk - fci) 

_ 2(fci-fc)fci |fci - fc|)|20H(fci - fc)5d(0? - 0?)/(fci, k, fci - fc) 

k 


/(A:) = 


47r^ 

/ic. 


poo 

r 

/ dki 

/ sin0id6(i 

Jo 

Jo 


(C2) 


with f{k, ki, ^ 2 ) = n{ki)n{k 2 ) — n{k 2 )n(k) — n{k)n{ki), 
01 = \^|A:—/cil/Zcb, and 02 = 'Ji{k+ki)/k\j. Finally, we 
can perform the integration by substitution with y = 0f 
and sin(0j) — Qj (j = 1, 2), obtaining Eq. (1351) . 


In summary, the convergence condition of the colli- 
sional integral is 4 < m < 5. In the Bogoliubov-wave 
turbulence, the power exponent is 9/2, so that this power 
law is local. 


Appendix D: Locality of the power law 


The locality condition of the power law is derived. In 
this derivation, we consider the convergence condition of 
collisional integral in Eq. ([35]). 

The power law is assumed to be n{k) oc with the 
positive power exponent m. Then, from Eq. (13511 . we find 
that the collisional integral may diverge at three points 
ki = 0, k, 00 . The point ki = k means the ^2 = 0, so that 
the convergence condition at ki = k is the same as that 
at fci = 0 because of the exchange symmetry (fci 0 / 82 ) 
of Eq. (E5)) . Thus, in the following, the convergence 
conditions at fci = 0, 00 are described. 

In the large fci region, the dangerous terms in the colli- 
sional integral are the second and third ones of Eq. (1351) , 
which is approximately expressed by 


2 



dki 


fcfn(fc) (ji{ki + k) 



n{ki — k 


(Dl) 


Thus, the convergence condition at fci = 00 is 4 < m. 

On the other hand, in the small fci region, the danger¬ 
ous terms appear in the first and third ones of Eq. ([M]) . 
These terms become 


k\n{ki)(^(k -I- fci) -I- n{k — fci) — 2n(fc)^ 


/ dki 

Jo 

f ^1 ,4 n ^^‘^n{k) 

oc [ dfcifc^-™, 

Jo 


(D2) 


which shows that the convergence condition at fci = 0 is 
m < 5. 
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